home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / see_r < prev    next >
Encoding:
Text File  |  1995-01-07  |  1.7 KB  |  73 lines
  1. %SEE    Pictures of a matrix and its (pseudo-) inverse.
  2. %       SEE(A) displays MESH(A), MESH(PINV(A)), SEMILOGY(SVD(A),'o'),
  3. %       and (if A is square) FV(A) in four subplot windows.
  4. %       SEE(A, 1) plots an approximation to the pseudospectrum in the
  5. %       third window instead of the singular values.
  6. %       SEE(A, -1) plots only the eigenvalues in the fourth window,
  7. %       which is much quicker than plotting the field of values.
  8. %       If A is complex, only real parts are used for the mesh plots.
  9. %       If A is sparse, just SPY(A) is shown.
  10.  
  11. require pinv fv ps rot90
  12.  
  13. see = function ( A, k )
  14. {
  15.   local (A, k)
  16.  
  17.   if (!exist (k)) { k = 0; }
  18.   m = A.nr; n = A.nc;
  19.   square = (m == n);
  20.  
  21.   B = pinv (A);
  22.   s = svd (A).sigma.';
  23.   zs = find (s == zeros (size (s)));
  24.   if (any( zs ))
  25.   {
  26.     % Remove zero singular values for semilogy plot.
  27.     s = s [complement (zs, 1:s.n)];
  28.   }
  29.  
  30.   if (length (ptmp = getplot ()) == 0)
  31.   {
  32.     pstart (2,2);
  33.   else if (ptmp.nplot < 4) {
  34.     error ("Create plot window with at least 4 sub-windows (pstart(2,2)).");
  35.   } }
  36.   
  37.   plaxis (); plimits (); plegend ();
  38.   Atmp = rot90 (A, -1);
  39.   plmesh (<< x=1:Atmp.nr; y=1:Atmp.nc; z=real (Atmp) >>);  # sub-plot 1
  40.  
  41.   plaxis (); plimits (); plegend ();
  42.   Btmp = rot90 (B, -1);
  43.   plmesh (<< x=1:Btmp.nr; y=1:Btmp.nc; z=real (Btmp) >>);  # sub-plot 2
  44.  
  45.   clear (Atmp, Btmp);
  46.  
  47.   if (k <= 0)
  48.   {
  49.     plgrid ("bcnst", "bcnstlv");
  50.     plstyle ("line-point");
  51.     plegend ();
  52.     plot (s);               # sub-plot 3
  53.  
  54.   else if (k == 1) {
  55.  
  56.     ps (A);                 # sub-plot 3
  57.  
  58.   } }
  59.  
  60.   if (square)
  61.   {
  62.     if (k == -1)
  63.     {
  64.       ps (A, 0);            # sub-plot 4
  65.     else
  66.       fv (A);               # sub-plot 4
  67.     }
  68.   else
  69.     printf ("Matrix not square.\n");
  70.   }
  71.   
  72. };
  73.